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^"V Abstract. The decrease of the barrier transparency in superconductor- 

insulator-superconductor (SIS) Josephson junctions leads to the devia- 
tions of the current-phase relation from the sinusoidal form. The sign of 
second harmonics is important for many applications, in particular in 
junctions with a more complex structure like SNINS or SFIFS, where N 
is a normal metal and F is a weak metallic ferromagnet. In our work we 
study the static magnetic flux distributions in long Josephson junctions 
£^ ■ taking into account the higher harmonics in the Fourier-decomposition of 

O | the Josephson current. Stability analysis is based on numerical solution 

^-» . of a spectral Sturm-Liouville problem formulated for each distribution. 

'_- ' In this approach the nullification of the minimal eigenvalue of this prob- 

Qh, lem indicates a bifurcation point in one of parameters. At each step of 

numerical continuation in parameters of the model, the corresponding 
nonlinear boundary problem is solved on the basis of the continuous 
analog of Newton's method. The solutions which do not exist in the tra- 
ditional model have been found. The influence of second harmonic on 
stability of magnetic flux distributions for main solutions is investigated. 

£^ ' Keywords: long Josephson junction, in-line geometry, Sturm-Liouville, 

O , double sine-Gordon, bifurcation, continuous analog of Newton's method, 

, ^*' l l fluxon, Numerov's finite-difference approximation 

^H ■ 

^ ! 1 Introduction 

\o 

t~ — ' Physical properties of magnetic flux in Josephson junctions (JJs) deserve the 

base of the modern superconducting electronics. Tunnel SIS JJs arc known to 
I/"") ■ be having the sinusoidal current phase relation. However, the decrease of the 

^D | barrier transparency in the SIS JJs leads the deviations of the current-phase 

relation from the sinusoidal form [TJ. We study the static magnetic flux distri- 
butions in the long JJs taking into account the second harmonic in the Fourier- 
decomposition of the Josephson current. The sign of the second harmonic de- 
pends on physical applications under considering. It is important, in particular, 
jJJ ■ in junctions like SNINS and SFIFS, where N is a normal metal and F is a weak 

Cd \ metallic ferromagnet [2]. Interesting properties of long Josephson junctions with 

an arbitrarily strong amplitude of second harmonic in current phase relation 
were considered in [3J. 
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Our purpose was to investigate an effect of the second harmonic accounting 
on the existence and stability magnetic flux distributions. Below, the numerical 
scheme and results of our stability analysis are demonstrated. 

2 Mathematical statement of the problem 

For a sufficiently wide class of JJ the superconducting Josephson current as a 
function of magnetic flux ip (phase difference of superconductors wave functions) 
can be represented as a sine series @|: 

oo 

Is = I c sin ip + y I m smrrup . (1) 



Using only first two terms of this expansion one can show [5] that the distribution 
of the magnitude (p(x) along x-axis of the junction in the static regime [4] satisfies 
the double sine-Gordon equation (2SG). 

— ip " + a± sin ip + a2 sin 2tp — 7 = , x € (— I; I) . (2) 

Here and below the prime means a derivative with respect to the coordinate x. 
The magnitude 7 is the external current, I is the semilength of the junction, 
ai and a,2 are parameters corresponding to I c and I2 in (JTJ respectively. They 
depend on the preparation technology of junctions |H6j . All the magnitudes are 
dimensionless. 

In the case of in-line geometry of the junction the boundary conditions for 
© have the form 

<p\±l) = K, (3) 

where h e is external magnetic field. 

From the mathematical viewpoint the transfer of the junction into dynamical 
regime [3] means |7l8j a stability loss (bifurcation) of all static solutions <p(x) of 
(J2J) , (El) a t the parameters 7 or h e variation. Our stability analysis of ip(x,p) was 
based on numerical solution of the corresponding Sturm-Liouville problem 

-i/>" + q(x)1) = \il), tp'(±l)=0 (4) 

with a potential q(x) = a,\ cosiy9 + 2a2 cos Ip. 

The minimal eigenvalue Ao(p) > corresponds the stable solution. In case 
Ao(p) < solution p(x,p) is unstable. The case Xq(j>) = indicates the bifurca- 
tion with respect to one of parameters p = (/, a±, 0,2, h e , 7). 



3 Numerical method 

Numerical solving of the boundary problem ©,© was performed on the basis of 
the Continuous analog of Newton's method [5] . At each Newtonian iteration the 
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corresponding linearized problem was solved using three-point Numerov's finite- 
difference approximation of the fourth order accuracy [§] . The discretization of 
the Sturm-Liouvillc problem (j4|) was realized with the help of standard second 
order finite-difference formulae. The calculation of the first several eigenvalues 
of the corresponding algebraic 3-diagonal problem was performed applying the 
standard subroutine from the package EISPACK. Details of numerical scheme 
are described in [TU] 
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Fig. 1. Change of Ao(7) for CS with increase of the coefficient 02 in the interval 
a 2 e [0; 0.7] at h e = 0, a t = 1, 21 = 10. 



4 Numerical results and conclusions 



Let us start with the trivial solutions of ([2])- hi the "traditional" case ai = 
two trivial solutions tp = and ip — n (below they are denoted by Mq and M v 
respectively) arc known at 7 = and h e = 0. Accounting of the second har- 
monic 02 sui2(p leads appearing two additional solutions ip = ± arccos(— a\j2ai) 
(denoted as M± ac ). The corresponding Ao as functions of 2SG-equation coeffi- 
cients have the form Ao[Mo] = a\ + 2a2, Ao[M T ] = —a\ + 2a2 and Ao[Af± ac ] = 
(af — 4a2)/2a2- The exponential stability of these constant solutions (CS) is 
determined by the signs of the parameters a\ and 02 and by its ratio a\/a% |10j . 

The dependencies of Ao on the external current 7 for CS at several positive 
values of 02 are demonstrated in Fig. [TJ Arising of the stable states M„. by the 
external current 7 at 02 > 0.5 is shown. 

When 02 < —0.5 the stable solution Mo disappears and other stable constant 
solutions M± ac arise. This transition is seen in Fig. [21 

Excepting CS, the 2SG equation is known to be supporting fluxon solutions. 
The fluxons play a significant role in the JJ physics. Different distributions of 
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Fig. 2. The same as on Fig. 1 but for a 2 € [-0.7; 0]. 
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Fig. 3. Internal magnetic field of the Fig. 4. The same as on Fig. 3 but 
fluxon & 1 at 7 = 0, h e — and 21 = 10 for decreasing a 2 < 0. 
when the parameter a 2 > increases. 



magnetic flux in J J are considered in the review [5]. At small external fields h e 
such distributions are fluxon $ l , antifluxon <P~ 1 and their bound states ^(p^ 1 
and <P~ l <P l . As external magnetic field h e is growing, more complicated stable 
fluxon and bound states appear: <? ±rl and ^ ±n ^T« ( n — \^ 2, 3, ...). 

Let us compare some basic physical characteristics of onc-fluxon solution 
# 1 in our model ©,([3]) with the traditional model (ai = 1, a 2 = 0). In Fig. 
[3] the deformation of the <p'(x) under influence of the parameter a 2 G [0; 1] 
is demonstrated. At a 2 = 0.5 the curve of internal magnetic field <p'(x) has a 
plateau in a neighborhood of the center x = 0. Further increase of the parameter 
a 2 leads to a formation of two maximums of the magnetic field. Thus, taking 
account of the coefficient a 2 leads to qualitative change of the form of fluxon 
distribution <2> . Such a change is not seen with a decrease in parameter a 2 when 
h e = (Fig. 2$. In the case of sufficiently large external magnetic field h e one 
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Fig. 5. Internal magnetic field of the Fig. 6. Full magnetic flux in depen- 

fluxon •P 1 at 7 = 0, h e = 2 and 21 = 10 dcncc on the parameter a 2 G [— 1; 1] 
when the parameter a 2 < decreases. at h e = 0, 7 = 0, 21 = 10 for <P l . 



can observe a similar qualitative deformation in the local minimums regions for 
a 2 < (see Fig.©. 

With change of the coefficient a 2 the number of fluxons [8] 



N(p) = —fv(x)dx, 



corresponding to the distribution & 1 is conserved i.e. dN/dci2 = 0. Here we have 
a value N^ 1 } = 1. 

When 02 is growing the full magnetic flux [8] A<p(p) = cp(l) — tp(—l) for this 
solutions tends to 2n, see Fig. [6] 

The value of the magnetic flux tp(x) in the middle of the interval does not 
change and ^(0) = 71*. 

In this paper we only focused on stability analysis of constant solutions and 
one-fluxon solutions in dependence on the 02 contribution. Investigation of an- 
other classes of solutions of 2GS-equation is the point of our further research. 
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